
############################################
############################################
# previous: 8-VAP_CDCut ####################
############################################

###### Summary ############################################################################################
# this script summarizes previous outputs of current and future habitat suitability into useful summary ###
# maps of patterns of change in species abundance, diveristy and turnover #################################
###########################################################################################################

############################################
# next: nothing, last script ###############
############################################

#################################################################################################################
#################################################################################################################
#################################################################################################################

#load packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(sp)
library(stringr)
library(gdistance)
library(Rcpp)
library(sf)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_occs")
SWDDir<-file.path(MyDir,"Maxent_SWDs/spp_swds")
BGDir<-file.path(MyDir,"Maxent_SWDs/backgrounds")
FullSDMDir<-paste(MyDir,"/Outputs1_Full",sep="")
OutDirCrossval<-paste(MyDir,"/Outputs2_Crossval",sep="")
OutDirFin<-file.path(MyDir,"Outputs3_Final")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")
IndSppOutDir<-paste(OutDirFin, "/Final_Individual_Spp", sep="")
HotSpotDir<-paste(OutDirFin, "/HotSpots",sep="")
MergeSppOutDir<-paste(OutDirFin, "/Final_Merge_Groups",sep="")

#set project name:
projectName<-"WHOSnakes"

#define selection criteria of interest for further analyses:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
critThresh<-c(1)                     #c(1,1,1,10,1,75)
futcats<-c("2050_Median","2090_Median","2050_Q10","2050_Q90","2090_Q10","2090_Q90")
#futcats<-c("2050_Median")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]
#spp<-DatSum$gen_sp_subtax
#spp<-DatSum$gen_sp_subtax[DatSum$PrimRegion=="Africa"]
#spp<-c("Bothrops_atrox","Bothrops_asper")


mapspp<-unique(DatSum$ModUnit[DatSum$Listing_comment!="not listed"])
mapspp<-paste(gsub(" ","_",mapspp))

bg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""), crs="+init=epsg:4326")
GLOBraster<-bg
values(GLOBraster)<-ifelse(values(GLOBraster)!=-Inf,0,-Inf) #raster showing all land as 0, ocean as NA
worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))

{
# #################################################################################
# #make temporary temp dir for rasters in this job and set from default to this dir
# dir.create(paste(FinalOutDir,"/",groupname,"_tmpCC",sep=""),showWarnings=FALSE)
# rasterOptions(tmpdir=paste(FinalOutDir,"/",groupname,"_tmpCC",sep=""))
# #################################################################################

############################################################################
#1. For each species, create rasters of increase/ decrease in range: #######
############################################################################

for (cat in futcats){
  for (sp in spp){

    if(

      (!file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_SuitIncrVsDecr.asc",sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_SuitIncrVsDecr.asc.gz",sep="")))
      |
      (!file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc.gz",sep="")))

       ){

    removeTmpFiles(h=0.05) #removes all temporary files older than 30min

  #calculate current species diversity
  if(file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc.gz", sep="")) &
     !file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc", sep="")) ){
  gunzip(filename=paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc.gz", sep=""), #unzip model
         destname=paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc", sep=""),
         remove=FALSE, overwrite=TRUE)
  }
  SDM1<-raster(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc", sep=""), crs="+init=epsg:4326")

  if(file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc.gz", sep="")) &
     !file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc", sep="")) ){
    gunzip(filename=paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc.gz", sep=""), #unzip model
         destname=paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc", sep=""),
         remove=FALSE, overwrite=TRUE)
  }
  SDM3<-raster(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc", sep=""), crs="+init=epsg:4326")

  myext<-extent(SDM1)




  #calculate range gain and loss areas for this species:
  print(paste(Sys.time(), "calculating potential range gain and loss for", sp, cat, sep=" "))

    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc.gz", sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc", sep="")) ){
     gunzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc.gz", sep=""), #unzip model
            destname=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc", sep=""),
            remove=FALSE, overwrite=TRUE)
    }
    SDM2<-raster(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc", sep=""), crs="+init=epsg:4326")
    speciesChange<- SDM2-SDM1 #both:1-1=0, future but not present: 1-0=1, present but not future: 0-1=-1
    writeRaster(speciesChange, filename=paste(IndSppOutDir,"/",sp, "_",cat,"_LossVsGain.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster

    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc", sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc.gz", sep=""))){
      gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc", sep=""),overwrite=TRUE, remove=TRUE)
    }

    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc", sep="")) &
       file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc.gz", sep=""))){
      file.remove(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc", sep=""))
    }

    if(file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc", sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_Current__CropThreshCDCut.asc.gz", sep=""))){
      gzip(filename=paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc", sep=""),overwrite=TRUE, remove=TRUE)
    }

    if(file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc", sep="")) &
       file.exists(paste(IndSppOutDir,"/",sp,"_Current__CropThreshCDCut.asc.gz", sep=""))){
      file.remove(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCutBin.asc", sep=""))
    }



    #calculate suitability increases and decreases this species:
    print(paste(Sys.time(), "calculating suitability changes for", sp, cat, sep=" "))

    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc.gz", sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc", sep="")) ){
     gunzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc.gz", sep=""), #unzip model
            destname=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc", sep=""),
            remove=FALSE, overwrite=TRUE)
    }
    SDM4<-raster(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc", sep=""), crs="+init=epsg:4326")
    speciesChange<- SDM4-SDM3 #both:1-1=0, future but not present: 1-0=1, present but not future: 0-1=-1
    writeRaster(speciesChange, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_SuitIncrVsDecr.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster



    #zip final files
    gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_SuitIncrVsDecr.asc",sep=""),overwrite=TRUE, remove=TRUE)
    gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep=""),overwrite=TRUE, remove=TRUE)




    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc", sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc.gz", sep=""))){
      gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc", sep=""),overwrite=TRUE, remove=TRUE)
    }

    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc", sep="")) &
       file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc.gz", sep=""))){
      file.remove(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc", sep=""))
    }

    if(file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc", sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc.gz", sep=""))){
            gzip(filename=paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc", sep=""),overwrite=TRUE, remove=TRUE)
    }

    if(file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc", sep="")) &
       file.exists(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc.gz", sep=""))){
      file.remove(paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc", sep=""))
    }

  }}
}





############################################################################
#2. Create raster of species gains, losses, turnover & richness change: ####
############################################################################
}

for (cat in futcats){
  
  TotalLosses<-GLOBraster
  TotalGains<-GLOBraster
  j<-spp[1]
  
  for (sp in spp[j:length(spp)]){
    
    n<-which(spp==sp)
    

    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc.gz",sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep="")) ){
      gunzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc.gz",sep=""), #unzip model
             destname=paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep=""),
             remove=FALSE, overwrite=TRUE)
    }
    
    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep=""))){
    speciesChange<-raster(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep=""), crs="+init=epsg:4326")
    speciesChange<-merge(speciesChange,GLOBraster)
    }else{
      speciesChange<-GLOBraster
      print(paste("loss vs. gain file does not exist for ",sp, ", continuing with previous version"))
    }
    
    print(paste("SpeciesChange:", sep=" "))
    print(speciesChange)
    
    print(paste(Sys.time(), "adding",sp,"to species gains & losses", sep=" "))
    
    ##############################
    #calculate losses
    if(file.exists(paste(HotSpotDir,"/",cat,"_TotalSpeciesLoss_prelim_",j-1,".asc",sep=""))){
      TotalLosses<-raster(paste(HotSpotDir,"/",cat,"_TotalSpeciesLoss_prelim_",j-1,".asc",sep="")) #read model output for species
    } else {
      TotalLosses<-TotalLosses
    }
    
    print(paste("TotalLosses before:", sep=" "))
    print(TotalLosses)
    
    #calculate cumulative species loss:
    print(paste(Sys.time(), "calculating cumulative range losses, adding ", sp, sep=" "))
    speciesLoss<-speciesChange
    values(speciesLoss)<-ifelse(values(speciesLoss)==(-1),1,0)
    
    print(paste("SpeciesLoss:", sep=" "))
    print(speciesLoss)
    
    TotalLosses<- TotalLosses+speciesLoss
    
    print(paste("TotalLosses after:", sep=" "))
    print(TotalLosses)
    
    writeRaster(TotalLosses, filename=paste(HotSpotDir,"/",cat,"_TotalSpeciesLoss_prelim_",j,".asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
    
    #remove previous backup file
    if(file.exists(paste(HotSpotDir,"/",cat,"_TotalSpeciesLoss_prelim_",j-2,".asc",sep=""))){
      file.remove(paste(HotSpotDir,"/",cat,"_TotalSpeciesLoss_prelim_",j-2,".asc",sep=""))
    }
    
    ###########################
    #calculate gains
    if(file.exists(paste(HotSpotDir,"/",cat,"_TotalSpeciesGain_prelim_",j-1,".asc",sep=""))){
      TotalGains<-raster(paste(HotSpotDir,"/",cat,"_TotalSpeciesGain_prelim_",j-1,".asc",sep="")) #read model output for species
    }else{
      TotalGains<-TotalGains
    }
    
    print(paste("TotalGains before:", sep=" "))
    print(TotalGains)
    
    #calculate cumulative species loss:
    print(paste(Sys.time(), "calculating cumulative range gains, adding ", sp, sep=" "))
    speciesGain<-speciesChange
    values(speciesGain)<-ifelse(values(speciesGain)==(1),1,0)
    
    print(paste("SpeciesGain:", sep=" "))
    print(speciesGain)
    
    TotalGains<- TotalGains+speciesGain
    
    print(paste("TotalGains after:", sep=" "))
    print(TotalGains)
    
    writeRaster(TotalGains, filename=paste(HotSpotDir,"/",cat,"_TotalSpeciesGain_prelim_",j,".asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
    
    #remove previous backup file
    if(file.exists(paste(HotSpotDir,"/",cat,"_TotalSpeciesGain_prelim_",j-2,".asc",sep=""))){
      file.remove(paste(HotSpotDir,"/",cat,"_TotalSpeciesGain_prelim_",j-2,".asc",sep=""))
    }
    
    ###################################
    #zip or remove original loss vs. gain file for this species:
    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep="")) &
       !file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc.gz",sep=""))){
      gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep=""),overwrite=TRUE, remove=TRUE)
    }
    
    if(file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep="")) &
       file.exists(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc.gz",sep=""))){
      file.remove(paste(IndSppOutDir,"/",sp,"_",cat,"_LossVsGain.asc",sep=""))
    }
    
    j<-j+1
  
    removeTmpFiles(h=0.1) #removes all temporary files older than 30min, then read hotspot map back in to continue process
    
    
    }
}

writeRaster(TotalLosses, filename=paste(HotSpotDir,"/",cat,"_TotalSpeciesLoss.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
writeRaster(TotalGains, filename=paste(HotSpotDir,"/",cat,"_TotalSpeciesGain.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster







print(paste(Sys.time(), "calculating overall species turnover and change in species richness", sep=" "))

for (cat in futcats){
  
  TotalLosses<- raster(paste(HotSpotDir,"/", cat, "_TotalSpeciesLoss.asc", sep=""), crs="+init=epsg:4326")
  TotalGains<-  raster(paste(HotSpotDir,"/", cat, "_TotalSpeciesGain.asc", sep=""), crs="+init=epsg:4326")
  CurrentSppRichness<-raster(paste(HotSpotDir,"/all_CurrentBin_HotSpot.asc",sep=""), crs="+init=epsg:4326")
  # FutureSppRichness<-raster(paste(HotSpotDir,"/all_",cat,"Bin_HotSpot.asc",sep=""), crs="+init=epsg:4326")
  
  #Spp Turnover   = (lost spp + gained spp ) / (total species (current plus new ones))
  SpeciesTurnover<- (TotalLosses+TotalGains) / (CurrentSppRichness + TotalGains)
  writeRaster(SpeciesTurnover, filename=paste(HotSpotDir,"/",cat,"_SppTurnover.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
  
  gzip(filename=paste(HotSpotDir,"/", cat, "_TotalSpeciesLoss.asc", sep=""),overwrite=TRUE, remove=TRUE)
  gzip(filename=paste(HotSpotDir,"/", cat, "_TotalSpeciesGain.asc", sep=""),overwrite=TRUE, remove=TRUE)
  gzip(filename=paste(HotSpotDir,"/",cat,"_SppTurnover.asc",sep=""),overwrite=TRUE, remove=TRUE)
  
  # #Change in richness = future richness - current richness
  # ChangeSppRichness  <- FutureSppRichness-CurrentSppRichness
  # writeRaster(ChangeSppRichness, filename=paste(HotSpotDir,"/",cat,"_SppRichnessChange.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
  
}





# ############################################################################
# #3. Create raster of individual & cumulative change in suitability:     ####
# ############################################################################
# 
# mergespp<-unique(DatSum$MergeSpecies)
# 
# for (cat in futcats){ #for each type of future: Q10,Q90, Median
#   
#   CurrentSppCumSuit<-raster(paste(HotSpotDir,"/all_CurrentSuit_HotSpot.asc",sep=""), crs="+init=epsg:4326")
#   FutureSppCumSuit<-raster(paste(HotSpotDir,"/all_",cat,"Suit_HotSpot.asc",sep=""), crs="+init=epsg:4326")
#   #Change in cumulative suitability = future cum suit - current cum suit
#   TotalSuitChange  <- FutureSppCumSuit-CurrentSppCumSuit
#   writeRaster(TotalSuitChange, filename=paste(HotSpotDir,"/",cat,"_TotalSuitabilityChange.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
#   
#   TotalSuitChange<-GLOBraster
#   
#   for (mergesp in mergespp){ #for each modelling unit (or should I do it for each merge group???)
#     
#     MaxSuitCurrent<-raster(paste(MergeSppOutDir,"/",mergesp,"_","CurrentSuit_Merged.asc",sep=""), crs="+init=epsg:4326") #read model output for species
#     MaxSuitFuture<-raster(paste(MergeSppOutDir,"/",mergesp,"_",cat,"Suit_Merged.asc",sep=""), crs="+init=epsg:4326") #read model output for species
#     MaxSuitChange<- MaxSuitFuture-MaxSuitCurrent
#     writeRaster(MaxSuitChange, filename=paste(MergeSppOutDir,"/",mergesp, "_", cat, "_SuitChange.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
#     
#     MaxBinCurrent<-raster(paste(MergeSppOutDir,"/",mergesp,"_","CurrentBin_Merged.asc",sep=""), crs="+init=epsg:4326") #read model output for species
#     MaxBiFuture<-raster(paste(MergeSppOutDir,"/",mergesp,"_",cat,"Bin_Merged.asc",sep=""), crs="+init=epsg:4326") #read model output for species
#     MaxBinChange<- MaxBinFuture-MaxBinCurrent
#     writeRaster(MaxBinChange, filename=paste(MergeSppOutDir,"/",mergesp, "_", cat, "_LossVsGain.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
#   }
#   
# }

# ######################################################
# #delete the temp dir for this script
# unlink(paste(FinalOutDir,"/",groupname,"_tmpCC",sep=""), recursive = TRUE) 
######################################################
######################################################################
######################################################################